11
2023-04-23
11 Image data
## ----data1, fig.cap=""--------------------------------------------------------
library("EBImage")
imagefile = system.file("images", "mosquito.png",
package = "MSMB")
mosq = readImage(imagefile)
## ----vis1, eval = FALSE-------------------------------------------------------
display(mosq)
## ----mosquito, fig.keep = 'high', fig.cap = "Mosquito discovered deceased in the suburbs of Decatur, Georgia (credit: CDC / Janice Haney Carr).", dev = "png", dpi = 300, fig.width=dim(mosq)[1]/300, fig.height=dim(mosq)[2]/300----
display(mosq, method = "raster")
text(x = 85, y = 800, label = "A mosquito",
adj = 0, col = "orange", cex = 1.5)
display(1-mosq, method = "raster")
## ----vis3a, eval = FALSE------------------------------------------------------
imagefile = system.file("images", "hiv.png",
package = "MSMB")
hivc = readImage(imagefile)
display(hivc)
## ----vis3b, eval = TRUE, echo = FALSE-----------------------------------------
imagefile = system.file("images", "hiv.png",
package = "MSMB")
hivc = readImage(imagefile)
## ---- hiv, eval = TRUE, echo = FALSE, fig.show = 'hold', fig.keep = 'high', fig.cap = "Scanning electron micrograph of HIV-1 virions budding from a cultured lymphocyte (credit: CDC / C. Goldsmith, P. Feorino, E.L. Palmer, W.R. McManus)."----
knitr::include_graphics(c('../images/hiv.png'), dpi = NA)
## ----image-oneminus, fig.keep = 'high', fig.cap = "Tiled display of four images of cell nuclei from the **[EBImage](https://bioconductor.org/packages/EBImage/)** package.", fig.margin = FALSE, dev = "png"----
nuc = readImage(system.file("images", "nuclei.tif",
package = "EBImage"))
nuc## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 510 510 4
## frames.total : 4
## frames.render: 4
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
## [2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
## [3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
## [4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
## [5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235dim(imageData(nuc))## [1] 510 510 4display(1 - nuc, method = "raster", all = TRUE)
display( nuc, method = "raster", all = TRUE)
## ----nucleitiledoneframe, dev = "png"-----------------------------------------
display(1 - nuc, method = "raster", frame = 2)
## ----how1---------------------------------------------------------------------
class(mosq)## [1] "Image"
## attr(,"package")
## [1] "EBImage"mosq[1:3,1:3]## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 3 3
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:3,1:3]
## [,1] [,2] [,3]
## [1,] 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000## ----how2, echo = FALSE, eval = FALSE-----------------------------------------
showClass("Image")## Class "Image" [package "EBImage"]
##
## Slots:
##
## Name: .Data colormode
## Class: array integer
##
## Extends:
## Class "array", from data part
## Class "structure", by class "array", distance 2
## Class "vector", by class "array", distance 3, with explicit coerce## ----dim----------------------------------------------------------------------
dim(mosq)## [1] 1400 952## ----mosqhist, fig.keep = 'high', fig.cap = "Histogram of the pixel intensities in `mosq`. Note that the range is between 0 and 1.", fig.width = .h$width, fig.height = .h$height----
hist(mosq)
## ----check, echo=FALSE--------------------------------------------------------
stopifnot(all(mosq>=0 & mosq<=1), isTRUE(all.equal(max(mosq), 1)), isTRUE(all.equal(min(mosq), 0)))
## ----how3---------------------------------------------------------------------
dim(imageData(mosq))## [1] 1400 952imageData(mosq)[1:3, 1:6]## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000 0.2039216 0.2000000 0.1960784## ----show1--------------------------------------------------------------------
mosq## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 1400 952
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:5,1:6]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [2,] 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784 0.1960784
## [3,] 0.1960784 0.1960784 0.2000000 0.2039216 0.2000000 0.1960784
## [4,] 0.1960784 0.1960784 0.2039216 0.2078431 0.2000000 0.1960784
## [5,] 0.1960784 0.2000000 0.2117647 0.2156863 0.2000000 0.1921569## ----show2--------------------------------------------------------------------
hivc## Image
## colorMode : Color
## storage.mode : double
## dim : 1400 930 3
## frames.total : 3
## frames.render: 1
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0## ----show3, echo=FALSE--------------------------------------------------------
stopifnot(colorMode(mosq)==Grayscale, colorMode(hivc)==Color, dim(nuc)[3]==4)
## ----how6, results = "hide"---------------------------------------------------
nuc## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 510 510 4
## frames.total : 4
## frames.render: 4
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
## [2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
## [3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
## [4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
## [5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235dim(imageData(nuc))## [1] 510 510 4## ----checkassertion, echo = FALSE---------------------------------------------
stopifnot(all(c(" frames.total : 4 ", " frames.render: 4 ") %in%
capture.output(EBImage:::showImage(nuc))))
## ----write1-------------------------------------------------------------------
writeImage(hivc, "hivc.jpeg", quality = 85)dim(hivc)## [1] 1400 930 3## ----objectsize---------------------------------------------------------------
object.size(hivc) |> format(units = "Mb")## [1] "29.8 Mb"object.size(hivc) |> format(units = "Kb")## [1] "30516.4 Kb"(object.size(hivc) / prod(dim(hivc))) |> format() |> paste("per pixel")## [1] "8 bytes per pixel"file.info("hivc.jpeg")$size## [1] 29490416 * 3 * 8 #three color, 16 Megapixel image## [1] 384## ---- mosqcrop, eval = TRUE, echo = FALSE, fig.show = 'hold', fig.keep = 'high', fig.cap = "Cropping: `mosqcrop`"----
knitr::include_graphics(c('../images/mosquito.png'), dpi = NA)
## ----manip1a------------------------------------------------------------------
mosqinv = normalize(-mosq)
## ----manip3a------------------------------------------------------------------
mosqcont = mosq * 3
mosqexp = mosq ^ (1/3)
display(mosqcont, method = "raster")
display(mosqexp, method = "raster")
display(mosq+0.4, method = "raster")
display(transpose(mosq), method = "raster")
mosq[100:438, 112:550]## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 339 439
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:5,1:6]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1882353 0.1882353 0.1921569 0.1882353 0.1882353 0.1882353
## [2,] 0.1921569 0.1921569 0.1921569 0.1960784 0.2000000 0.1921569
## [3,] 0.1960784 0.1960784 0.1960784 0.2078431 0.2078431 0.1960784
## [4,] 0.1960784 0.1882353 0.1843137 0.1960784 0.2039216 0.2039216
## [5,] 0.1960784 0.1803922 0.1764706 0.1882353 0.2000000 0.2078431## ----manip4a------------------------------------------------------------------
mosqcrop = mosq[100:438, 112:550]
mosqthresh = mosq > 0.5
mosqtransp = transpose(mosq)
display(mosqcrop, method = "raster")
display(mosqthresh, method = "raster")
display(mosqtransp, method = "raster")
## ----checkassertionont, echo = FALSE------------------------------------------
stopifnot(identical(t(mosq), transpose(mosq)))## ----spattrans1---------------------------------------------------------------
mosqrot = EBImage::rotate(mosq, angle = 30)
mosqshift = translate(mosq, v = c(40, 70))
mosqflip = flip(mosq)
mosqflop = flop(mosq)
display(mosqrot, method = "raster")
display(mosqshift, method = "raster")
display(mosqflip, method = "raster")
display(mosqflop, method = "raster")
## ----MSMB, results="hide"-----------------------------------------------------
imagefiles = system.file("images", c("image-DAPI.tif",
"image-FITC.tif", "image-Cy3.tif"), package="MSMB")
cells = readImage(imagefiles)
display(cells)## Only the first frame of the image stack is displayed.
## To display all frames use 'all = TRUE'.
## ----checkdim, echo=FALSE-----------------------------------------------------
stopifnot(dim(cells)[3]==3)
## ----range--------------------------------------------------------------------
apply(cells, 3, range)## image-DAPI image-FITC image-Cy3
## [1,] 0.001586938 0.002899214 0.001663233
## [2,] 0.031204700 0.062485695 0.055710689## ----fixrange-----------------------------------------------------------------
cells[,,1] = 32 * cells[,,1]
cells[,,2:3] = 16 * cells[,,2:3]
apply(cells, 3, range)## image-DAPI image-FITC image-Cy3
## [1,] 0.05078202 0.04638743 0.02661173
## [2,] 0.99855039 0.99977111 0.89137102getFrame(cells, 1)## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 340 490
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:5,1:6]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.05517662 0.05420005 0.05566491 0.05371176 0.05517662 0.05322347
## [2,] 0.05566491 0.05517662 0.05420005 0.05420005 0.05566491 0.05468833
## [3,] 0.05371176 0.05615320 0.05517662 0.05517662 0.05371176 0.05322347
## [4,] 0.05517662 0.05468833 0.05420005 0.05566491 0.05420005 0.05273518
## [5,] 0.05566491 0.05468833 0.05664149 0.05371176 0.05517662 0.05371176## ----defw---------------------------------------------------------------------
w = makeBrush(size = 51, shape = "gaussian", sigma = 7)
nucSmooth = filter2(getFrame(cells, 1), w)
dim(w)## [1] 51 51## ----image-filter2, fig.keep = 'high', fig.cap = "The middle row of the weight matrix, `w[` (ref:image-filter2-1) `, ]`", fig.width = 3, fig.height = 2.6, dev = "png"----
library("tibble")
library("ggplot2")
tibble(w = w[(nrow(w)+1)/2, ]) |>
ggplot(aes(y = w, x = seq(along = w))) + geom_point()
## ----smooth-------------------------------------------------------------------
cellsSmooth = Image(dim = dim(cells))
sigma = c(1, 3, 3)
for(i in seq_along(sigma))
cellsSmooth[,,i] = filter2( cells[,,i],
filter = makeBrush(size = 51, shape = "gaussian",
sigma = sigma[i]) )## ----illuminationartifact1----------------------------------------------------
py = seq(-1, +1, length.out = dim(cellsSmooth)[1])
px = seq(-1, +1, length.out = dim(cellsSmooth)[2])
illuminationGradient = Image(
outer(py, px, function(x, y) exp(-(x^2+y^2))))
nucBadlyIlluminated = cellsSmooth[,,1] * illuminationGradient## ----illuminationartifact2----------------------------------------------------
disc = makeBrush(21, "disc")
disc = disc / sum(disc)
localBackground = filter2(nucBadlyIlluminated, disc)
offset = 0.02
nucBadThresh = (nucBadlyIlluminated - localBackground > offset)## ----adathresh----------------------------------------------------------------
nucThresh =
(cellsSmooth[,,1] - filter2(cellsSmooth[,,1], disc) > offset)
## ----morphopen1---------------------------------------------------------------
nucOpened = EBImage::opening(nucThresh,
kern = makeBrush(5, shape = "disc"))## ----imageProcessing14--------------------------------------------------------
nucSeed = bwlabel(nucOpened)
table(nucSeed)## nucSeed
## 0 1 2 3 4 5 6 7 8 9 10
## 155408 511 330 120 468 222 121 125 159 116 520
## 11 12 13 14 15 16 17 18 19 20 21
## 115 184 179 116 183 187 303 226 164 309 194
## 22 23 24 25 26 27 28 29 30 31 32
## 148 345 287 203 379 371 208 222 320 443 409
## 33 34 35 36 37 38 39 40 41 42 43
## 493 256 169 225 376 214 228 341 269 119 315## ----imageProcessing17, eval = FALSE------------------------------------------
## display(colorLabels(nucSeed))
## ----imageProcessing15a-------------------------------------------------------
nucMask = cellsSmooth[,,1] - filter2(cellsSmooth[,,1], disc) > 0
## ----imageProcessing15b-------------------------------------------------------
nucMask = fillHull(nucMask)
## ----imageProcessing16--------------------------------------------------------
nuclei = propagate(cellsSmooth[,,1], nucSeed, mask = nucMask)## ----voronoiExample-----------------------------------------------------------
zeros = Image(dim = dim(nuclei))
voronoiExamp = propagate(seeds = nuclei, x = zeros, lambda = 100)
voronoiPaint = paintObjects(voronoiExamp, 1 - nucOpened)
## ----voronoiEx----------------------------------------------------------------
head(table(voronoiExamp))## voronoiExamp
## 1 2 3 4 5 6
## 5645 4735 370 5964 3333 1377ind = which(voronoiExamp == 13, arr.ind = TRUE)
head(ind, 3)## row col
## [1,] 112 100
## [2,] 113 100
## [3,] 114 100## ----histcellbody-1, fig.width = .h$width, fig.height = .h$height, fig.keep="high", eval = TRUE, echo = FALSE, fig.cap = "Histogram of the actin channel in `cellsSmooth`, after taking the logarithm."----
hist(log(cellsSmooth[,,3]) )
## ----histcellbody-2, fig.width = .h$width, fig.height = .h$height, fig.keep="high", eval = TRUE, echo = FALSE, fig.cap = "Zoom into Figure \\@ref(fig:histcellbody-1)."----
hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)
## ----histcellbody, eval = FALSE, echo = TRUE----------------------------------
## hist(log(cellsSmooth[,,3]) )
## hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)
## ----checkhistcellsSmooth, echo=FALSE-----------------------------------------
stopifnot(mean(cellsSmooth[,,3]>=exp(-3.6) & cellsSmooth[,,3]<=exp(-3.1)) > 0.68)
## ----musigmaEstimator, message=FALSE------------------------------------------
library("genefilter")
bgPars = function(x) {
x = log(x)
loc = half.range.mode( x )
left = (x - loc)[ x < loc ]
wid = sqrt( mean(left^2) )
c(loc = loc, wid = wid, thr = loc + 6*wid)
}
cellBg = apply(cellsSmooth, MARGIN = 3, FUN = bgPars)
cellBg## [,1] [,2] [,3]
## loc -2.90176965 -2.94427499 -3.52191681
## wid 0.00635322 0.01121337 0.01528207
## thr -2.86365033 -2.87699477 -3.43022437## ----histcellbody-3, fig.keep = 'high', fig.cap = "As in Figure \\@ref(fig:histcellbody-2), but with `loc` and `thr` shown by vertical lines.", fig.width = .h$width, fig.height = .h$height----
hist(log(cellsSmooth[,,3]), xlim = -c(3.6, 3.1), breaks = 300)
abline(v = cellBg[c("loc", "thr"), 3], col = c("brown", "red"))
## ----cytoplasmMask------------------------------------------------------------
cytoplasmMask = (cellsSmooth[,,2] > exp(cellBg["thr", 2])) |
nuclei | (cellsSmooth[,,3] > exp(cellBg["thr", 3]))
## ----imageProcessing22--------------------------------------------------------
cellbodies = propagate(x = cellsSmooth[,,3], seeds = nuclei,
lambda = 1.0e-2, mask = cytoplasmMask)## ----imageProcessing25--------------------------------------------------------
cellsColor = rgbImage(red = cells[,,3],
green = cells[,,2],
blue = cells[,,1])
nucSegOnNuc = paintObjects(nuclei, tgt = toRGB(cells[,,1]),
col = "#ffff00")
nucSegOnAll = paintObjects(nuclei, tgt = cellsColor,
col = "#ffff00")
cellSegOnAll = paintObjects(cellbodies, tgt = nucSegOnAll,
col = "#ff0080")## ----baserfeats---------------------------------------------------------------
meanNucInt = tapply(cells[,,1], nuclei, mean)
meanActIntInNuc = tapply(cells[,,3], nuclei, mean)
meanActIntInCell = tapply(cells[,,3], cellbodies, mean)
## ----pairsint, fig.keep = 'high', fig.cap = "Pairwise scatterplots of per-cell intensity descriptors. \\label{pairsint}", fig.margin = FALSE, fig.width = 5.3, fig.height = 5.3----
library("GGally")## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(tibble(meanNucInt, meanActIntInNuc, meanActIntInCell))
## ----imageProcessing27--------------------------------------------------------
F1 = computeFeatures(nuclei, cells[,,1], xname = "nuc",
refnames = "nuc")
F2 = computeFeatures(cellbodies, cells[,,2], xname = "cell",
refnames = "tub")
F3 = computeFeatures(cellbodies, cells[,,3], xname = "cell",
refnames = "act")
dim(F1)## [1] 43 89## ----showF1-------------------------------------------------------------------
F1[1:3, 1:5]## nuc.0.m.cx nuc.0.m.cy nuc.0.m.majoraxis nuc.0.m.eccentricity nuc.0.m.theta
## 1 119.5523 17.46895 44.86819 0.8372059 -1.314789
## 2 143.4511 15.83709 26.15009 0.6627672 -1.213444
## 3 336.5401 11.48175 18.97424 0.8564444 1.470913## ----readlymphnodedata--------------------------------------------------------
library("readr")##
## Attaching package: 'readr'## The following object is masked from 'package:genefilter':
##
## speclibrary("dplyr")##
## Attaching package: 'dplyr'## The following object is masked from 'package:EBImage':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncellclasses = c("T_cells", "Tumor", "DCs", "other_cells")
brcalymphnode = lapply(cellclasses, function(k) {
read_csv(file.path("..", "data",
sprintf("99_4525D-%s.txt", k))) |>
transmute(x = globalX,
y = globalY,
class = k)
}) |> bind_rows() |> mutate(class = factor(class))## Rows: 103681 Columns: 5## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 27822 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 878 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 77081 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): filename
## dbl (4): locX, locY, globalX, globalY
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.brcalymphnode## # A tibble: 209,462 × 3
## x y class
## <dbl> <dbl> <fct>
## 1 6355 10382 T_cells
## 2 6356 10850 T_cells
## 3 6357 11070 T_cells
## 4 6357 11082 T_cells
## 5 6358 10600 T_cells
## 6 6361 10301 T_cells
## 7 6369 10309 T_cells
## 8 6374 10395 T_cells
## 9 6377 10448 T_cells
## 10 6379 10279 T_cells
## # ℹ 209,452 more rowstable(brcalymphnode$class)##
## DCs other_cells T_cells Tumor
## 878 77081 103681 27822## ----checkcellnumber, echo = FALSE--------------------------------------------
tabtab = table(brcalymphnode$class)
within = function(x, a, b) (x>a & x<b)
stopifnot(all(within(tabtab[c("T_cells", "Tumor", "DCs")], c(100000, 27000, 800), c(110000, 28000, 1000))))
## ----brcalntcells, fig.keep = 'high', fig.cap = "Scatterplot of the $x$ and $y$ positions of the T- and tumor cells in `brcalymphnode`. The locations were obtained by a segmentation algorithm from a high resolution version of Figure \\@ref(fig:stainedlymphnode). Some rectangular areas in the T-cells plot are suspiciously empty, this could be because the corresponding image tiles within the overall composite image went missing, or were not analyzed.", fig.margin = FALSE, dev = "png", dpi = 300, fig.width=9, fig.height=4.5, pointsize=24----
ggplot(filter(brcalymphnode, class %in% c("T_cells", "Tumor")),
aes(x = x, y = y, col = class)) + geom_point(shape = ".") +
facet_grid( . ~ class) + guides(col = "none")
## ----spatstat1, results = "hide"----------------------------------------------
library("spatstat")## Loading required package: spatstat.data##
## Attaching package: 'spatstat.data'## The following object is masked _by_ '.GlobalEnv':
##
## cells## Loading required package: spatstat.geom## spatstat.geom 3.0-3##
## Attaching package: 'spatstat.geom'## The following object is masked from 'package:genefilter':
##
## area## The following objects are masked from 'package:EBImage':
##
## affine, closing, distmap, opening, rotate## Loading required package: spatstat.random## spatstat.random 3.0-1## Loading required package: spatstat.core## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## Loading required package: rpart## spatstat.core 2.4-0## Loading required package: spatstat.linnet## spatstat.linnet 2.3-2##
## spatstat 2.3-3 (nickname: 'That's not important right now')
## For an introduction to spatstat, type 'beginner'## ----spatstat2----------------------------------------------------------------
ln = with(brcalymphnode,
ppp(x = x, y = y, marks = class, xrange = range(x), yrange = range(y)))
ln## Marked planar point pattern: 209462 points
## Multitype, with levels = DCs, other_cells, T_cells, Tumor
## window: rectangle = [3839, 17276] x [6713, 23006] units## ----checkclassln, echo = FALSE-----------------------------------------------
stopifnot(identical(class(ln), "ppp"))## ----spatstat3----------------------------------------------------------------
library("spatstat")
ln = with(brcalymphnode, ppp(x = x, y = y, marks = class,
xrange = range(x), yrange = range(y)))
ln## Marked planar point pattern: 209462 points
## Multitype, with levels = DCs, other_cells, T_cells, Tumor
## window: rectangle = [3839, 17276] x [6713, 23006] units## ----densityppz1--------------------------------------------------------------
d = density(subset(ln, marks == "Tumor"), edge=TRUE, diggle=TRUE)
plot(d)
## ----densityppp1, fig.keep = 'high', fig.cap = "Intensity estimate for the cells marked `Tumor` in `ppp`. The support of the estimate is the polygon that we specified earlier on (Figure \\@ref(fig:convhull)).", fig.width = 3.75, fig.height = 3.5, echo = FALSE----
par(mai = c(0, 0, 0.2, 0.7))
plot(d)
## ----densityppz0--------------------------------------------------------------
d0 = density(subset(ln, marks == "Tumor"), edge = FALSE)
plot(d0)
## ----densityppp0, fig.keep = 'high', fig.cap = "As Figure \\@ref(fig:densityppp1), but without edge correction \\label{densityppp0}", fig.width = 3.75, fig.height = 3.5, echo = FALSE----
par(mai = c(0, 0, 0.2, 0.7))
plot(d0)
## ----relrisk-calc-------------------------------------------------------------
rr = relrisk(ln, sigma = 250)
## ----relrisk, fig.keep = 'high', fig.cap = "Estimates of the spatially varying probability of each of the cell claases, conditional on there being cells.", fig.margin = FALSE, fig.width=7, fig.height=7----
plot(rr)
## ----checkConditional, echo=FALSE---------------------------------------------
m = rr[[1]]$v
for(i in 2:length(rr)) m = m + rr[[i]]$v
#stopifnot(all(is.na(m) | abs(m-1)<1e-6))
## ----Gestshow-----------------------------------------------------------------
gln = Gest(ln)
gln## Function value object (class 'fv')
## for the function r -> G(r)
## .....................................................................
## Math.label Description
## r r distance argument r
## theo G[pois](r) theoretical Poisson G(r)
## han hat(G)[han](r) Hanisch estimate of G(r)
## rs hat(G)[bord](r) border corrected estimate of G(r)
## km hat(G)[km](r) Kaplan-Meier estimate of G(r)
## hazard hat(h)[km](r) Kaplan-Meier estimate of hazard function h(r)
## theohaz h[pois](r) theoretical Poisson hazard function h(r)
## .....................................................................
## Default plot formula: .~r
## where "." stands for 'km', 'rs', 'han', 'theo'
## Recommended range of argument r: [0, 21.033]
## Available range of argument r: [0, 61.89]## ----Gest, fig.keep = 'high', fig.cap = "Estimates of $G$, using three different edge effect corrections --which here happen to essentially lie on top of each other-- and the theoretical distribution for a homogenous Poisson process. \\label{Gest}", fig.width = 4.25, fig.height = 4.25, results="hide"----
library("RColorBrewer")
plot(gln, xlim = c(0, 10), lty = 1, col = brewer.pal(4, "Set1"))
## ----Linhom-------------------------------------------------------------------
Lln = Linhom(subset(ln, marks == "T_cells"))## number of data points exceeds 1000 - computing border correction estimate onlyLln## Function value object (class 'fv')
## for the function r -> L[inhom](r)
## ................................................................................
## Math.label
## r r
## theo L[pois](r)
## border {hat(L)[inhom]^{bord}}(r)
## bord.modif {hat(L)[inhom]^{bordm}}(r)
## Description
## r distance argument r
## theo theoretical Poisson L[inhom](r)
## border border-corrected estimate of L[inhom](r)
## bord.modif modified border-corrected estimate of L[inhom](r)
## ................................................................................
## Default plot formula: .~.x
## where "." stands for 'bord.modif', 'border', 'theo'
## Recommended range of argument r: [0, 819.84]
## Available range of argument r: [0, 819.84]## ----Images-Lln, fig.keep = 'high', fig.cap = "Estimate of $L_{\\scriptsize \\mbox{inhom}}$, Equations \\@ref(eq:kinhom) and \\@ref(eq:Lest), of the T cell pattern. \\label{Images-Lln}", fig.width = 4, fig.height = 5.9, results = "hide"----
plot(Lln, lty = 1, col = brewer.pal(3, "Set1"))
## ----pcfdo--------------------------------------------------------------------
pcfln = pcf(Kinhom(subset(ln, marks == "T_cells")))## number of data points exceeds 1000 - computing border correction estimate only## ----Images-pcf, fig.keep = 'high', fig.show = "hold", fig.cap = "Estimate of the pair correlation function, Equation \\@ref(eq:pcf), of the T cell pattern. \\label{Images-pcf}", fig.width=5, fig.height=5, results="hide"----
plot(pcfln, lty = 1)
plot(pcfln, lty = 1, xlim = c(0, 10))
## ----samplingpcf, fig.keep = 'high', fig.cap = "Answer to Question \\@ref(ques:Image-ques-samplingpcf): as in the bottom panel of Figure \\@ref(fig:Images-pcf), but with denser sampling.", fig.width=5, fig.height=5----
pcfln2 = pcf(Kinhom(subset(ln, marks == "T_cells"), r = seq(0, 10, by = 0.2)))## number of data points exceeds 1000 - computing border correction estimate onlyplot(pcfln2, lty = 1)